Superconductivity and charge order of confined Fermi systems 
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The low-temperature properties of the two-dimensional attractive Hubbard model are strongly 
influenced by the fermion density. Away from half-filling, there is a finite-temperature transition to a 
phase with s-wave pairing order. However, T c is suppressed to zero at half-filling, where long-range 
charge-density-wave order also appears, degenerate with superconductivity. This paper presents 
Determinant Quantum Monte Carlo simulations of the attractive Hubbard model in the presence 
of a confining potential Vtrap which makes the fermion density p inhomogeneous across the lattice. 
Pair correlations are shown to be large at low temperatures in regions of the trapped system with 
incommensurate filling, and to exhibit a minimum as the local density p(i) passes through one 
fermion per site. In this ring of p(i) — 1, charge order is enhanced. A comparison is made between 
treating Vtrap within the local-density approximation (LDA) and in an ab initio manner. It is argued 
that certain sharp features of the LDA result at integer filling do not survive the proximity of doped 
sites. The critical temperature of confined systems of fixed characteristic density is estimated. 

PACS numbers: 67.85.-d, 05.30.Fk, 05.10.Ln 



INTRODUCTION 

Studies of the interplay of spatial inhomogeneity 
and superconductivity have a long history. A sem- 
inal early result was Anderson's realization^ that al- 
though the breaking of translation invariance by dis- 
order renders momentum no longer a good quantum 
number, pairing still occurs between appropriately cho- 
sen (time reversed) states. Numerical studies within 
the Bogoliubov-de Gennes approximation; 2 ^ Quantum 
Monte Carlo (QMC)^, and other approaches have quan- 
tified the magnitude of disorder which superconductivity 
can withstand 5 -. In these studies, and the granular su- 
perconducting materials they model 6 -, regions of pairing 
order coexist with normal, or insulating, phases. Super- 
conductivity can be destroyed by various mechanisms, 
including phase fluctuations between the order param- 
eter on different islands where Cooper pairs of bosons 
existji or breaking of the Cooper pairs themselves^ The 
dominant mechanism determines the appropriate mod- 
eling, e.g. a description within the disordered boson^ or 
fermion Hubbard Hamiltonians, or "phase-only" descrip- 
tions with the XY model and its variants.— 

Recently, experiments on ultra-cold atoms have pro- 
vided a rather different realization of inhomogeneity 
in the form of a smoothly varying confining potential 
which produces a system with a radial density profile, 
maximal at the trap center and falling off to zero at 
the periphery^ Much attention has focussed on repul- 
sively interacting bosons and fermions j 12 ' 13 In this case, 
a Mott insulator may coexist with supcrfluid or nor- 
mal phases. For fermions, the Mott insulator also ex- 
hibits antifcrromagnctic correlations at very low temper- 



atures. At present, experimentally accessible tempera- 
tures for fermionic systems are such that a degenerate 
Fermi gas has been observed^, along with signatures 
of the Mott phasei 15 ' 16 The ultimate objective is insight 
into the ground state physics of the repulsive Hubbard 
model (RHM), and, in particular, the fundamental issue 
of d-wave superconducting order and its interplay with 
antifcrromagnetismjii 

This goal for repulsive fermions awaits the attainment 
of lower experimental temperatures. In the interim, it 
is useful to perform careful studies of attractive systems. 
This case is not only of interest in its own right, but also 
QMC simulations can often attain lower temperatures for 
attractive models, and thus can track experiments closer 
to transitions into ordered phases. 

The focus of the present paper is the description of the 
behavior of attractively interacting fermions in a two- 
dimensional confining potential. Some of the issues are 
similar to the repulsive case, in particular the coexistence 
of phases as the density varies across the trapJ^ However, 
the attractive case has several important distinctions, 
specifically the existence of known finite-temperature 
phase transitions in two dimensions. In addition, in the 
repulsive case there is a broad range of chemical poten- 
tials p which fall within the "Mott gap" and for which 
the fermion density p = 1. That is, the compressibil- 
ity k = dp/dp = at p = 1. For the confined system, 
this implies an extended region of commensurate density, 
spatial sites which have a value of the local confining po- 
tential which falls within the Mott gap. In the attractive 
case, the compressibility is finite (n ^ 0) at commensu- 
rate density. As a consequence, the region of half-filling 
is a truly one-dimensional ring as opposed to an annulus 
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of finite thickness. 

A key result of this work is that the unique features of 
charge density wave physics at the single value of chem- 
ical potential which gives commensurate filling do not 
survive coupling to neighbors of incommensurate density. 
Thus the correlations which appear in a homogeneous 
system with commensurate filling are never achieved in a 
trap; the Local-Density Approximation (LDA), in which 
the behavior of each site in a confining potential is as- 
sumed to be that of a homogeneous system with global 
density matching the local filling, breaks down at that 
point. 

This paper is organized as follows: In the next sec- 
tion we describe the specific Hamiltonian, the attrac- 
tive Hubbard model ( AHM) and aspects of the computa- 
tional methodology, Determinant Quantum Monte Carlo 
(DQMC), which will be used. Results are then presented 
within the LDA as well as from direct simulations of con- 
fined systems, and the two approaches are compared. 
Next, we present a finite-size extrapolation using data 
from systems of different sizes at constant characteristic 
density, and estimate the critical temperature of the con- 
fined AHM. A concluding section summarizes the results 
and indicates some remaining open questions. 

Studies of the AHM with inhomogeneity have been 
performed with Variational Monte Carlo,— Bogoliubov 
de-GenneS ) 20 ' 21 and Gutzwillcr approaches ! 22 ' 23 Of par- 
ticular relevance here is work within dynamical mean 
field theory (DMFT) and a two-site impurity solver—, 
which suggested that the half-filled physics is stabilized 
by a confining potential, and that an extended supersolid 
phase of commensurate density exists in a trap. 



MODELS AND COMPUTATIONAL APPROACH 

The attractive Hubbard Hamiltonian, in the presence 
of a confining potential, reads 

(*i> i 

~Y.it 1 ~ Krap \i\ 2 } (n n + nil ) . (1) 
i 

Here c ia (c io .) are creation (destruction) operators at spa- 
tial site i for two different species of fermions a. We 
choose the center of the trap to be at a plaquette cen- 
ter and set the origin there, so that the coordinates i x 
and i y take half-integer values. In the condensed-matter 
context, a = ±i is the electron spin. For cold atoms, a 
labels two hypcrfine states. We will consider the case of 
square lattices of linear size L. The hopping parameter t 
can be tuned by changing the optical lattice depth,—; in 
the following, t = 1 is chosen to set the scale of energy. 
The sum Y)uj\ is over all near- neighbor pairs of sites, 
and J^i is over a U sites. The on-site attraction \U\ can 
be tuned through the application of a magnetic field via 



a Fcshbach resonance. The chemical potential /i is set to 
get the desired number of particles TV. Finally, Vtrap is 
the trap curvature which determines the strength of the 
confining potential. 

In DQMGA the partition function Z = Tver 1311 is 
written as a path integral by discrctizing the inverse tem- 
perature /3 into M intervals of size At = f3/M. The 
Trotter approximation^ e~ ArH ~ e~ e~ v isolates 
the quartic terms (involving the interaction U) in H, and 
a discrete Hubbard-Stratonovich field^ decouples e~ ArV 
so that only quadratic terms in the fermion operators ap- 
pear. When the trace over fermion operators is done, Z 
is expressed as a sum over the different field configura- 
tions with a weight which is the product of two determi- 
nants (one for each value of a) of matrices with dimension 
L 2 x L 2 given by the number of lattice sites. In the case 
of attractive U, because the two species couple to the 
Hubbard-Stratonovich field with the same sign, the two 
determinants are identical and there is no sign problem. 2 - 
This allows us to study confined systems down to ar- 
bitrarily low temperatures, unlike the repulsive model 
where the largest j3 accessible is /3 ~ 3-4 for confined 
systems with U = 4-8^ 

The observables which will be the focus of this paper 
are the s-wave pairing and charge density wave (CDW) 
correlation functions, 

Cpair(i, j) = (A i+J -At), ^ 
c c hargo(i,i) = (ni+jUi) - {n i+j } (rii) , 

where A] = cLc], creates a pair of fermions on site i 
and rii = + counts the fermions on site i. Notice 
that these depend on i and not just on the separation j. 
We also define the associated structure factors 

"Scdw = ^ c chargc(*7 j)- 

ij 

In addition, we study the local quantities 
p(i) = (n it ) + (nn), 

(ij) 

*(i) = !E( c L^+4 c -)> (4) 

3 

D(i) = (nifnii) - (THf){nii); 

the density p(i) has already been used; t(i) is the kinetic 
energy associated with the bonds of site i; and D(i) is 
the double occupancy with the trivial density dependence 
subtracted. 

Before proceeding with the confined case, it is useful 
to review the properties of the translationally invariant 
case, V trap = 0. In two dimensions, it is known that 
the half-filled attractive Hubbard Hamiltonian has com- 
bined long-range CDW and s-wave pairing order in its 
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ground state, and is unordered at nonzero temperature*^ 
When doped, the symmetry between charge and pairing 
is broken, and a Kosterlitz-Thouless (KT) transition to 
a quasi-long-range ordered superfluid phase occurs at fi- 
nite temperature. Many numerical and analytical stud- 
ies have been performed.— The transition temperature 
T c rises rapidly as p is made non-zero and reaches a 
maximum value of T c ~ 1 /w for a wide range of fill- 
ings 0.5 < p < 0.9j22 The effect of inhomogencitics in the 
interaction strength has also been explored^ 

Consideration of the "asymmetric" particle-hole (PH) 
transformation helps clarify these assertions. On a bipar- 
tite lattice, when the down-spin operations in the AHM, 
Eq. Q] are mapped with 

the kinetic energy is unchanged, but the interaction term 
changes sign, so that the AHM maps onto the RHM (the 
phase factor (— l) 1 is understood to take the values ±1 
on alternating sublattices). This PH symmetry provides 
a simple argument that the two-dimensional half-filled 
AHM can have long-range order (LRO) only at T = 0, 
like the RHM. 

QMC simulations have shown that the ground state of 
the half-filled, two-dimensional uniform RHM is magnet- 
ically ordered^ PH symmetry then implies that CDW 
and pair order occur simultaneously in the T = half- 
filled AHM. To see this, note that the z component of spin 
jii-f — nil m the RHM maps onto the charge + in 
the AHM, so that magnetic LRO in the z direction of the 
RHM corresponds to CDW order of the AHM. Similarly, 
magnetic order in the xy plane maps onto ,s-wavc pairing 
order. The degeneracy of the z and xy magnetic order 
in the repulsive model implies that CDW and pair order 
occur simultaneously in the half-filled attractive case. 

A final consequence of PH symmetry is the explana- 
tion of the occurrence of pairing order (and the absence 
of CDW order) at finite temperature in the doped AHM. 
When doped, p is non-zero. Under the PH transforma- 
tion Eq.[5j the chemical potential term p(ni^+riii) in the 
AHM becomes a Zeeman field B{n^ — n^) in the RHM 
(B has the same numerical value as fx). Because the or- 
der in the RHM is antiferromagnetic, a uniform field in 
the z direction makes it energetically favorable for spins 
to lie in the xy plane, since then they can tilt out of the 
plane and pick up field energy without costing as much 
exchange energy. This lowering of symmetry from three 
to two components makes possible a finite temperature 
Kosterlitz-Thouless transition in two dimensions. The 
xy magnetic order which exists in the RHM then maps 
to s-wave pair order in the AHM. 

CORRELATIONS AND THE LOCAL DENSITY 
APPROXIMATION 

We begin by showing the density profile in Fig.[TJ along 
with the local kinetic energy and the double occupancy. 




FIG. 1. (Color online) Fermion density p, kinetic energy t 
and double occupancy D as a function of distance |i| from 
the trap center, with TV = 564 particles, U = 6, Kra P = 
0.0097 and p. = 0.8 at (3 = 9. The solid line shows the LDA 
result; in this case, near perfect agreement is found between 
the true trapped system and the LDA. Note that, in contrast 
to the RHM, there is no plateau at half-filling (marked by the 
vertical line), which corresponds to the absence of a Mott gap 
in the homogeneous model. Error bars are smaller than the 
symbols. 



Results are given both from the LDA and from a trapped 
30 x 30 system with N = 564.1 ± 0.4 particles. These 
two approaches yield results in very good agreement for 
p(i), t(i) and D(i). An important point is the absence 
of a density plateau at p = 1 , in accordance with results 
in the LDA, and also with a particle-hole symmetry ar- 
gument which identifies the compressibility n = dp/ dp 
of the AHM with the uniform magnetic susceptibility 
X = dM/dB of the RHM, which is known to be nonzero 
at zero external field; thus, as noted above, the AHM 
has finite compressibility at integer filling, and there is 
no Mott plateau at half-filling. 

This true one-dimensionality of the p(i) = 1 ring makes 
the formation of long range CDW order in the AHM 
much less robust than the antiferromagnetic order which 
can occur on the quasi-two-dimensional Mott annulus of 
integer filling that occurs in the RHM. 

A related difference to the RHM is seen in the kinetic 
energy, which shows a maximum at half-filling in Fig. [TJ 
where in the RHM the Mott phase would lead to localiza- 
tion and a minimum of the kinetic energy. This behav- 
ior is best understood by applying the asymmetric PH 
transformation. The corresponding RHM is uniformly 
half-filled and subject to a perpendicular Zeeman field Bi 
that varies radially and goes through zero at those sites 
i that had p(i) = 1 in the original, attractive model. 
Away from the B — region, the system gets increas- 
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FIG. 2. (Color online) Near-neighbor [c pa ir (i, (J))] alm next- 
near neighbor [c pa ir (i, ( } ))] s-wave pair correlations as a func- 
tion of distance |£| from the trap center, with iV = 564 parti- 
cles, (7 = 6, Krap = 0.0097 and p = 0.8 at p = 9. The solid 
lines show the LDA result; in this case, a striking failure of 
the LDA is seen around half-filling (marked by the vertical 
line). Due to CDW/s-wave degeneracy at half- filling in the 
homogeneous model, the LDA predicts superfluidity only suf- 
ficiently far from that point, while in the trapped system the 
superfiuid phase may penetrate the half-filled ring. Error bars 
have been suppressed to avoid clutter, but the spread of the 
data points gives an indication of the uncertainties, which re- 
sult from large statistical fluctuations observed at low T. (See 
text.) 



ingly spin-polarized and thereby makes Pauli exclusion 
more effective in hindering fermion mobility. One there- 
fore expects a maximum in both \t(i)\ and |-D(i)| when 
the site i belongs to the Bi = region. As these quan- 
tities are unchanged by the PH transformation, this last 
statement translates verbatim to the AHM. 

Next, Fig. [2] shows the near- neighbor c pa ; r (i, ( J )) and 
next-near neighbor c pa i r (i, ( \ )) s-wave pairing correla- 
tors both in the LDA and the 30 x 30 system. In the 
LDA, both functions dip at r = ju/^trap; where the 
local density p{%) = 1, as do the corresponding farthest- 
neighbor correlators. Fig. [3] shows the density correla- 
tors charge (*, ( o)) an d |c c harge («, ( } )) which in the LDA 
peak as the system crosses through commensurate filling. 

The dip (peak) in the pairing (CDW) correlation func- 
tions observed in the LDA may be understood from the 
CDW-pairing degeneracy that exists precisely at half- 
filling, and the corresponding suppression of the finite-T 
pairing order that exists away from half-filling. 

While the LDA compares favorably with the ab ini- 
tio calculation across most of the lattice, the dip (peak) 
in the s-wave pairing (CDW) when the p(i) = 1 ring is 
crossed is conspicuously absent in the true trapped sys- 
tem. 
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FIG. 3. (Color online) Near-neighbor [c pa i r (i, (J))] an< i next- 
near neighbor [c pa ir (i, ( i )), with the sign inverted for clar- 
ity] CDW correlations as a function of distance \i\ from the 
trap center, with N = 564 particles, [7 = 6, Vtrap = 0.0097 
and p = 0.8 at f3 — 9. The solid lines show the LDA re- 
sult. As with the pair correlation, the LDA fails around 
half-filling (marked by the vertical line), where it predicts 
enhanced CDW correlations. Error bars are smaller than the 
symbols. 



It is useful to compare this behavior with the RHM, 
where the physics of commensurate filling (p = 1) can 
be inferred correctly, for the most part, from the LDA 
because of the presence of an annulus of finite thickness 
which "protects" the Mott region. By contrast, in the 
AHM, there is no such protection; the half-filled ring is 
truly one-dimensional, and the physics of commensurate 
filling is essentially absent in the trap. 

Once again, a deeper understanding can be reached 
by applying the PH transformation: the B = region 
of the corresponding RHM divides, and is coupled to, 
regions where the spins are tilted out of the xy plane in 
opposite directions. Spins in the -6 = region can then 
lower the system's energy by aligning with neighboring 
spins on the xy plane and therefore breaking the local 
5(7(2) symmetry characteristic of the Hubbard model in 
the absence of an external field. In the original language 
of the AHM this implies that, at half-filling, one should 
expect a reduction of the CDW correlation, an increase 
in the pairing correlations and a breaking of the CDW- 
pairing degeneracy. 

For optical-lattice experiments that aim to emulate the 
Hubbard model, our findings indicate that the physics of 
the half-filled AHM will be inaccessible in any experimen- 
tal setup that leads to a confining potential (fi— Vtrap |*| 2 ) 
as in Eq. [T] It has recently been suggested^ to simulate 
the AHM on an optical lattice, and then utilize the PH 
transformation Eq. [S] to draw conclusions about the re- 
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pulsive case. However, this proposal appears to be chal- 
lenging, since the effective absence of the half-filled case 
in the confined AHM means that the physics of the RHM 
at zero Zeeman field will be inaccessible. 

Both problems outlined above could be solved by the 
putative "off-diagonal" confinement (ODC)24, if realized, 
since particle-hole symmetry indicates that the inhomo- 
gencous lattice can be made uniformly half-filled under 
ODC. On the other hand, in the conventional "diagonal 
confinement" , observing a finite-temperature transition 
to a superconducting phase in confined systems becomes 
much more likely when the CDW region only occupies 
such a limited spatial region, as discussed in detail in the 
following section. 



FINITE-SIZE EXTRAPOLATION 

We now turn to the interesting issue of the finite-T 
phase transition in the confined AHM. In general we 
may ask, when a trapping potential is added to a model 
that undergoes, e.g. a KT transition in the homogeneous 
case, how is the nature of that phase transition altered 
by the trap? In the case of the classical xy-model, it 
has been showr*2£ that the KT transition of the homoge- 
neous model is preserved in many respects in the trapped 
case. In this section, we present a finite-size extrapola- 
tion (FSE) to address the same question for the AHM. 

True phase transitions of course can occur only in the 
thermodynamic limit of infinite system size. In a trans- 
lationally invariant system, the correct way to perform 
this limit is familiar and almost trivial: the global den- 
sity p = N/L 2 is kept constant. In the presence of a trap, 
the correct generalization is to keep the the "characteris- 
tic density" p := N/£ 2 constant^ Here I = y^/Krap is 
the natural length scale in the problem, formed by com- 
bining the kinetic energy t and the trap curvature Vtrap- 
In the finite-size analysis described below, we have used 
systems of constant p. 

Of course, since only finite system sizes are accessible 
numerically, one must invoke a procedure to infer be- 
havior in the thermodynamic limit from finite lattices. 
Here, we follow the finite-size scaling (FSS) procedure 
discussed in [32| . The generalization of FSS in the pres- 
ence of the trapping potential has been called trap-size 
scaling (TSS)2£; we note that TSS may be expressed as 
keeping p constant and proceeding with the FSS analysis 
as usual. 

In this approach, the pair structure factor Ps of Eq. [3] 
is obtained for different lattice sizes and temperatures 
(plotted in Fig. [4j. We see that at high temperatures 
{P < 4), when the correlation length is short, Ps is in- 
dependent of the system size. At lower temperatures 
((3 ~ 4), the curves begin to separate as the correla- 
tion length becomes large compared to the system sizes. 
Finally, for very low temperatures Ps approaches a con- 
stant depending on the system size. This behavior is 
expected for KT and second order transitions where ob- 
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FIG. 4. (Color online) The s-wave structure factor Ps for 
systems of different linear size L with U = 6 as a function 
of inverse temperature /3. At low temperatures (/3 > 7), the 
curves flatten off to values determined by the system size, 
indicating a divergent correlation length. 

servables stop evolving with temperature when the cor- 
relation length £ exceeds L. 

Unfortunately the full FSS procedure of [Hf is de- 
feated in the case of the confined AHM because of ex- 
cessive statistical fluctuations in the DQMC estimator 
for c pa j r . These large fluctuations are manifest in, for 
example, the very large error bars at in two of the near- 
neighbor pair correlators of Fig. [2] However, a compar- 
ison of results on lattice sizes 12, 14, . . . , 24 at tempera- 
tures 1 < j3 < 16 provides evidence that the KT transi- 
tion is indeed preserved in the trap. The approach hinges 
on the critical scaling of the structure factor near the 
phase transition; the expected KT form is 

P S ~ L 2 -^ (6) 

with r](T c ) = 1 /a. The critical exponent 77 for a KT tran- 
sition in a homogeneous system is known to vary with 
temperature between ri(T = 0) = and rj(T c ) = 
In the trapped system, this issue is complicated by the 
varying filling. Arguing within the LDA, since the filling 
varies in the lattice, so does T c ; but 77 must be a func- 
tion of T/T c rather than T itself, therefore r\ should vary 
across the system along with the filling. We must there- 
fore ask whether an effective exponent 77 exists such that 
Eq. [6] holds for the system as a whole. 

To investigate the behavior of r\ in the trapped system, 
we plot log Ps as a function of log L for several low tem- 
peratures, j3 > 10, in Fig. \5\ A straight line with slope 
2 — i) in the doubly logarithmic plot indicates that the 
scaling relation Eq. [5] holds for the whole system with 
an effective exponent 77. We have used non-linear least- 
squares fits to estimate r) at temperatures /3 > 6.6, using 
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FIG. 5. (Color online) The s-wave structure factor Ps as 
a function of linear lattice size L for U — 6 at various low 
temperatures (top to bottom: /3 = 1 /t — 16, 15, . . . , 10) on a 
doubly logarithmic scale. A vertical shift has been applied to 
separate the curves from each other. A linear dependence on 
the doubly logarithmic scale (i.e., the structure factor varies 
as a power of the system size) points to a divergent correlation 
length. 

the data shown in Fig. 01 Within the statistical uncer- 
tainties, we find good agreement with Eq. [6] some exam- 
ples of the fits are shown in Fig. [5] 

As a function of temperature, the behavior of fj is con- 
sistent with the behavior expected from the KT tran- 
sition in the homogeneous case, dropping off to zero as 
the temperature is decreased. Between (3 = 8-10, we find 
large fluctuations of fj, which may be interpreted as a sig- 
nature of the phase transition. Thus, we may estimate 
the critical temperature of the system as T c ~ 0.1-0.125. 

Our results for fj are consistent with recent work on 
the classical xy-mode^, where it was found that fj for 
the trapped system at a certain temperature takes the 
same value as 77 in the homogeneous system at the same 



temperature. 



CONCLUSIONS 

In this paper, the attractive Hubbard model (AHM) 
in a harmonic confining potential was studied, especially 
with regards to superfluidity at low temperatures. Re- 
sults from the local-density approximation (LDA) were 
compared to calculations within the true trapped system. 

While the LDA is valid for local observablcs in most 
cases, we find qualitatively wrong predictions for the s- 
wave pairing and CDW correlation functions around the 
ring of half-filling, where the LDA predicts a dip and a 
peak, respectively, which are absent when the trapped 
system is treated ab initio. This is linked to the relation- 
ship between the pairing and CDW correlations at half- 
filling; and to the density profile in the trap, where the 
Mott plateau exhibited by the repulsive Hubbard model 
is absent. Consequently, the physics of the half-filled 
case is not represented anywhere in the trap. On the 
other hand, that the CDW correlations are suppressed 
will make it easier to observe a transition to a superfluid 
phase at finite temperature. 

A finite-size extrapolation, where systems of different 
sizes but at the same characteristic density p are com- 
pared, provides evidence that the Kostcrlitz-Thouless 
transition to a superfluid phase of the homogeneous 
AHM persists in the trap, with a critical temperature 
of T c ~ 0.1-0.125. 
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